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We argue for the existence of a liquid ground state in a class of square lattice models of orbitally 
degenerate insulators. Starting with the SU(4) symmetric Kugel-Khomskii model, we utilize a Ma- 
jorana Fermion representation of spin-orbital operators to access novel phases. Variational wave 
functions of candidate liquid phases are thus obtained, whose properties are evaluated using vari- 
ational Monte Carlo. These states are disordered, and are found to have excellent energetics and 
ground state overlap (> 40%) when compared with exact diagonalization on 16 site clusters. We 
conclude that these are spin-orbital liquid ground states with emergent nodal fermions and Z2 gauge 
fields. Connections to spin 3/2 cold atom systems and properties in the absence of SU(4) symmetry 
are briefly discussed. 



I. INTRODUCTION 



In correlated insulators, the degrees of freedom that 
remain at low energies are spin and orbital degeneracy. 
At low temperatures one usually obtains an ordered state 
described essentially by a classical variable, the Landau 
order parameter. Ground states that are not described by 
the Landau framework are expected to possess strikingly 
new properties. While they are known to occur in one 
dimensional systems, an important question is whether 
they arise in bulk two and three dimensional materials. 
Theoretical studies have largely focused on quantum spin 
systems. While model Hamiltonians for spin liquids exist, 
one needs special conditions like strong frustration to en- 
sure that the spins do not order. Otherwise, even for spin 
1/2 quantum fluctuations are not typically strong enough 
to destroy order. On the other hand, orbital degener- 
acy in insulators can enhance quantum fluctuations jL2i2i 
destroying order and possibly lead to a spin-orbital liq- 
uid state. An experimental illustration is provided by 
the insulating spinels MnSc2S4 and FeSc2S4. The for- 
mer, a pure S=5/2 system, magnetically orders below 
2Kelvin. The latter, a spin S=2 system, is identical in 
most respects except that it involves a twofold orbital 
degeneracy. In contrast to the spin system, it is found to 
remain spin disordered down to the lowest temperatures 
of 30mKelvin. 4 

Here, we theoretically study a simple spin 1/2 square 
lattice model with the minimal twofold orbital degen- 
eracy. Such a spin-orbital Hilbert space is realized for 
e.g. with the d 3z 2_ r 2, d x 2_ y 2 orbitals of Ag 2+ , Cu 2+ or 
(low spin) Ni 3+ in an octahedral environment. We focus 
on a model that captures the effect of enhanced quan- 
tum fluctuations from orbital degeneracy, and argue for 
a liquid ground state in this case. We give theoretical ar- 
guments for the stability of this phase as well as well as a 
numerical Monte Carlo study of a variational wave func- 
tion. The latter is found to have extremely good ener- 
getics for our Hamiltonian, and allows us to characterize 
this phase beyond the simple fact that it is disordered. 
The low energy collective excitations of the liquid state 



are captured by an emergent Z2 gauge field, coupled to 
Dirac like fermionic excitations with fractional spin and 
orbital quantum numbers. The ground state wavefunc- 
tion is strongly entangled, and can be thought of as a 
product of three Slater determinant wave- functions. 

Realistic spin-orbital Hamiltonians tend to be rather 
complicated with several exchange couplings that are 
strongly direction dependent. Moreover, a linear cou- 
pling between the orbital degrees of freedom and the 
Jahn- Teller phonons can quench coherent orbital dynam- 
ics. However, for sufficiently strong exchange interac- 
tions, the coupling to phonons can be ignored, and the 
orbitals can be taken to be quantum degrees of freedom. 
Since we are interested here in the general effects of en- 
hanced quantum fluctuations from orbital degeneracy, we 
follow 3 and others in considering a model that treats all 
four states on a site symmetrically, i.e. the SU(4) sym- 
metric spin orbital model. This will allow for compar- 
isons and is a useful starting point. We show later that 
our essential conclusions are unchanged on perturbing 
away from this high symmetry point. 

The high symmetry SU(4) model may, in fact, have a 
direct physical realization. We point out that a model 
of e g orbitals on certain high symmetry lattices, like the 
diamond lattice, can have spin orbital Hamiltonian that 
are nearly SU(4) symmetric. A different setting for this 
physics has been opened up by the recent experimen- 
tal developments on the trapping and cooling of alkaline 
earth atoms. Confining these to the sites of an opti- 
cal lattice leads to SU(N) symmetric magnetic models. 
The nuclear spin here provides the N flavors, and, given 
the weak dependence of scattering lengths on nuclear 
spin, leads to SU(N) symmetric exchange interactions, 
which, for fermionic atoms, will be antiferromagnetic£ 
Fermionic alkaline-earth like atoms of 173 Yb have been 
cooled to quantum degeneracy^ while the Mott state of 
the bosonic 174 Yb has been recently realizedi 7 - A differ- 
ent realization may be provided by spin 3/2 cold atoms, 
such as 132 Cs confined to the sites of an optical lattice. 
At unit filling, one has four states per site, and although 
the physical symmetry is only that of spin rotations, the 
small difference in scattering lengths imply only a weak 
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breaking of SU(4) symmetry. It was pointed out in£ that 
even including these differences only breaks the symme- 
try down to SO(5)xZ2 in the low energy limit. 



II. THE MODEL: 

We study the SU(4) symmetric Kugel-Khomskii 
model&i£ on the square lattice: 

<ij> 

where s are the spin- 1/2 Pauli matrices and r are the 
Pauli matrices acting on the two degenerate orbital 
states. We consider the antiferromagnetic(AFM) case 
(J > 0) and set J = I hereafter. The high symmetry of 
this model implies that the three spin operators s, three 
orbital operators f and nine spin orbital operators a a r b 
all appear with equal weight. These are fifteen genera- 
tors of SU(4). It should be noted that symmetry does 
not uniquely define a model, one also needs to specify 
the representation of the symmetry group appearing at 
each site. Here the fundamental representation appears 
and an SU(4) singlet can only be formed between four 
sites. 

The model ([1]) was numerically studied m£ using exact 
diagonalization (the model suffers from a sign problem in 
the spin-orbital basis) on system sizes of up to 4x4. The 
ground state (E = —17.4) is an SU(4) singlets at zero 
wavevector. Simple trial states, such as with spin-orbital 
order , or a box singlet state, have much higher average 
energy (E = and E — —12 respectively) pointing to the 
importance of quantum fluctuations. This motivates the 
study of spin orbital liquids as candidate ground states. 

The idea of resonating valence bonds ; 11 ' 12 provides an 
intuitive picture of the quantum spin liquid. A more 
formal approach that is easier to generalize is the slave 
particle formalism, where the spin is decomposed into 
'partons', e.g. s r = a' a ar& a a'r> wnere the a are the 
Pauli matrices, and (a| r , aj r ) creates a boson (Schwinger 
boson) or fermion with spin (t, 4) at site r, and the con- 
straint aJ. r a CTr = 1 is imposed at every site. One then 
makes a mean field decomposition to obtain a quadratic 
Hamiltonian, and the constraint is then imposed by pro- 
jecting the wave function to obtain a variational ground 
state. The state so obtained is a candidate spin-liquid 
wave function. This procedure can be generalized in a 
straightforward way to the spin-orbital Hilbert space at 
hand, by introducing a four component a a and the con- 
straint above at every sitei 13 i 14 The physical spin and or- 
bital operators are again bilinears of a a . However, there 
are some drawbacks to this straightforward generaliza- 
tion. The bosonic parton representation cannot treat the 
SU(4) symmetric point, while the fermionic parton theo- 
ries necessarily lead to Fermi surface states which can be 
hard to stabilize as ground states. 



Below, we will sidestep these difficulties by showing 
that this problem admits a third, physically distinct, par- 
ton representation in terms of Majorana Fermions. This 
representation, which has not previously been applied to 
two dimensional systems, offers us many advantages. Be- 
sides being more economical (in terms of expanding the 
Hilbert space in the minimal fashion), it leads to liquid 
states with a Z2 gauge group, whose low energy physics is 
well understood and know to exist as stable phases. We 
emphasize here that the projected wavefunction obtained 
from this Majorana parton representation of the SU(4) 
model is distinct from the 'Schwinger' fermion represen- 
tation, involving four fermionic a a operators. 



III. MAJORANA PARTON FORMULATION: 

We first point out the group isomorphism 
SU(4)=SO(6). The 15 generators of the latter are 
dimension six antisymmetric real matrices L^ v , where 
A = 1 ... 15 and /1, v = I ... 6. An operator repre- 
sentation of this algebra is obtained by introducing 
six Majorana fermions (xi, xe) which satisfy the 
anticommutation relations {x^j Xv\ — 2(5 M „. The 
operators A = jL^XftXvt where summation over 
repeated indices is assumed, reproduce the commutation 
relations for SO (6) generators. The Majorana Fermions 
transform as SO (6) vectors. 

We now use the group isomorphism to obtain a repre- 
sentation of the spin-orbital operators, in terms of Majo- 
rana fermions. It is helpful to write the set of Majorana 
fermions as a pair of three component vectors (9 r , ff r ) 
where, e.g. 9 r — (9i r , Oir, &3r) and we have introduced 
site indices r. The spin and orbital operators can then 
be written in the compact form: 

s r = - -ff r x ff r , T r = --d r x9 r (2) 

and s^t" = —irj^O^r, which automatically obey the ex- 
pected algebra. Note, the sign of the Majorana fermion 
operators can be changed without affecting the physical 
operators. This Z2 redundancy is connected to the fact 
that the Hilbert space is now enlarged - since each Ma- 
jorana fermion corresponds to \/2 degrees of freedom, we 
have (\/2) 6 = 8 states whereas there are only 4 physical 
states per site. The excess states can be removed by im- 
plementing a Z2 constraint at each site: first define the 
operator v r commutes with the physical operators (which 
are fermion bilinears) and is idempotent t£ = 1. 



Hence implementing the constraint in ([3]), restricts us 
to the physical Hilbert space. This operator generates 
the Z 2 gauge transformation on the Majorana fermions 
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The model in ([I]) can be written in these variables as: 



H 



E 

<jk> 



(1/8) ( 



(4) 



When supplemented by the constraint ([3]), this is an ex- 
act rewriting of the model. Here, the SO(6) symmetry of 
the model is explicit. The quartic nature of the Hamilto- 
nian requires an approximation. We begin with a mean 
field treatment and use it to generate variational states 
in which the constraint is treated exactly. 

In the context of models with only spin 1/2 (and no 
orbital degrees of freedom) we note that a representa- 
tion utilizing three Majorana fermions per site, where 
the spin operator s r is given by an expression identi- 
cal to that in Eqn. [2J has been studied*^ However, we 
point out that this is distinct from our current formalism 
since a single site constraint cannot be applied to gen- 
erate the physical Hilbert space. Nevertheless, this pro- 
vides an alternate parton approach - for example, with 
an even number of sites, one can view half of the spins 
as 'orbital pseudo-spins', then the spin problem will be 
artificially converted to a spin-orbital problem, (although 
without SU(4) symmetry). Then the gauge fixing, and 
construction of complex fermions and variational wave 
functions can be proceeded in the same fashion as in this 
paper. But this artificial discrimination of spin and 'or- 
bital pseudo-spin' usually will superficially break lattice 
symmetry. Another way to construct a Hilbert space for 
the spin 1/2 only model, is to introduce a fourth Majo- 
rana fermion on every site, and the set of four fermions 
satisfies a product constraint as in Eqn.O 1 ^ This has the 
benefit of being formulated with a unique single site con- 
straint. However, unfortunately it turns out that these 
four fermions are just the real and imaginary parts of the 
two Schwinger fermion operators, and do not lead to a 
new representation. The constraint is the familiar one of 
requiring single occupancy of the Schwinger fermions. 

Mean Field Theory and Gutzwiller Projection: 
With real mean field parameters \jk (— ~Xkj)i we have: 



H 



MF 



E 

<jk> 



[l - (i Xjfe /4) • ffk + 6j ■ Ok) + X%/8 



J5) 

In self consistent mean field theory \jk — HiVj ' Vk + 
Sj ■ 9k)) mf- For convenience we combine the 6 Majorana 



fermions into 3 complex fermions: 



-t = 



(l/2)fa 



iO ar ) which are more intuitive although the SO(6) sym- 
metry is no longer explicit. The constraint then is: 
2 a =i 4r c ar = OR 2 (while the odd values of the site 
occupation are forbidden). Writing c r = (c\ r , C2 r , c^ r ), 
we have i(ffj ■ ffk + 9j ■ 9k) = 2i(c] • Ck — 4 • Cj) i.e. the 
mean field theory simply involves fermions hopping with 
pure imaginary amplitudes. Such a band structure is au- 
tomatically particle-hole symmetric, which leads to half- 
filled bands for each of the c ra fermions. Note, despite 
the imaginary hoppings the mean field ansatz is time re- 
versal symmetric if the hopping is bipartite. The mean 



field wave function is simply a product of three identical 
Slater determinants. While the specific Slater determi- 
nant depends on the mean field ansatz, we make a few 
general observations below. If we consider a system with 
47V sites, required to obtain an SU(4) singlet state, each 
Slater determinant $ is a function of 27V particle co- 
ordinates, corresponding to half filling. Gutzwiller pro- 
jecting the mean field state into the constrained Hilbert 
space, yields a physical spin-orbital wave-function. In the 
fermion representation, a site can either have no fermions 
(denoted by |0)), or two fermions, in which case there are 
three states, \X) = c|4|0), \Y) = c|c||0>, \Z) = 44l )- 
These are related to the spin orbital basis states via: 

\a* = T l,r 2 =±1) = (|0>±i|X»A/2 
\a z = ±l,T* = ±l) = (\Y)±i\Z))/V2 

Given a configuration specified by the locations of the 
\X), \Y) and \Z) states (at sites {xi},{yj} and {z m } re- 
spectively, where Xi, yj, z m are 37V distinct positions), 
the spin-orbital wave function assigns an amplitude 
^ [{xi}, {yi}, {zi}} to it. Note, the locations of the |0) 
states are automatically specified. For an SU(4) singlet 
we need equal numbers, TV, of the four types of sites, so 
{x^ = {x\, . . . , xn} etc. After the Gutzwiller projection 
we obtain: 



{z m }} ■ $[{z m }, { Xi }] ■ $[{ Xi }, {y 3 }} 



(6) 



Thus the projected spin-orbital wave function is a prod- 
uct of three Slater determinants with a lot of entangle- 
ment. We now apply this formalism to specific models. 



IV. ONE DIMENSIONAL CHAIN 

This SU(4) symmetric nearest neighbor model in ID 
is very well understood and serves as a benchmark for 
our technique. The only symmetric mean field ansatz 
is a uniform Xr,r+i = X, leading to a dispersion e(k) = 
xsin(fc). We construct the resulting projected wave func- 
tion for L site chain with L = 8, 16 . . . 128 and antiperi- 
odic boundary conditions, and evaluate its properties us- 
ing variational Monte Carlo. The energy per site from 
the projected wave functions extrapolated to the ther- 
modynamic limit is —0.8233, not far from the exact 

result, 17 1 - (1/2) ^^ dx = -0.8251. The lead- 
ing term in the asymptotic spin correlation function is 
cos(7rr/2)/r 15 , consistent with theoretical and numeri- 
cal predictions i 18 ' 19 i 20 Note, these desireable properties 
of the wavefunction only arise after projection. (see 
TABLE U. 

Interestingly, the n/2 wavevector of the dominant cor- 
relations do not correspond to a natural wavevector of the 
mean field dispersion, and arises entirely from projection. 
In contrast, this wavevector is easier to understand on 
projecting a quarter filled band, which arises in the stan- 
dard fermionic representation of spin-orbital operators 
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TABLE I: Results of the projected wave function for L- 
site chain with antiperiodic boundary condition. The sec- 
ond row shows energy per site for the projected wave func- 
tion. The third row shows the L -1 ' 5 scaling of s z -correlation 
functions [s z (a;) is s z at position x]. 



L 


8 16 32 48 64 


(H)/L 


-0.8642 -0.8332 -0.8256 -0.8242 -0.8237 


(s z (0)s z (L/2))L 1 ' 5 


1.66 1.82 1.90 1.92 1.96 


L 


80 96 112 128 


(H)/L 


-0.8235 -0.8234 -0.8233 -0.8233 


{s%Q)s'(L/2))L 1 - 5 


1.98 1.94 1.97 1.97 



i v 
A i A 



A 



FIG. 1: The 7r-flux ansatz on the square lattice. An arrow 
from site j to k means Xjfc > 0. Dashed lines enclose the 
doubled unit cell, with two sites u and v. 



a a . Remarkably, one can show that the projected wave 
functions arising from this representation and the Ma- 
jorana fermion representation discussed above, are iden- 
tical in one dimension. We stress that this is a special 
feature of one dimension, and in higher dimensions, the 
two will lead to physically distinct states. Details of the 
proof can be found in Appendix [A"l 



V. SQUARE LATTICE: 

The mean field states on the square lattice can be dis- 
tinguished by the gauge invariant flux through the ele- 
mentary plaquettes e.g. XjkXkiXimXmj for the plaque- 
tte jk£m. Translation and Time reversal symmetry dic- 
tates that this flux must be uniform and can be either 
or 7T. This leads to two distinct mean field states the 
uniform and 7r flux state ansatz. The uniform states 
ansatz is Xr.r+st = X?,r+y = X, where r = (x, y) is 
the position of lattice sites. The mean field dispersion 
is e(k) = xi sm (kx) + s'm(k y )] for all the three flavors, and 
has a square Fermi surface. However, the uniform ansatz 
state has higher energy than the 7r-flux ansatz both in 
mean field theory and after projection, so we focus on 
the 7r-flux ansatz. 

The 7r-flux ansatz is (— l) v Xr,r+x — Xr,r+y — Xi as 
shown in FIG. [TJ The unit cell in mean field theory is 
doubled, with u and v sublattices as shown. 

After a Fourier transform, the mean field Hamiltonian 



© is: 

Hmf = (2 + X 2 /4)L 2 

+xY(c f , c t e ) ( si 

^ L. / \ a,k,u a,/c,u / I gj 

a,k \ 



sin k x 
sin k v 



sin k 




where the sum over k is over the L 2 /2 (for L x L lattice) 
fc-points in the reduced(0 < k x < 2ir, < k y < ir) Bril- 
louin zone(BZ), and a = {1,2,3}. The above result can 
be further diagonalized by a Bogoliubov transformation 
and produce the two branches of the mean field disper- 
sion e±(fc) = ±x^/ sin 2 k x + sin 2 k y . This dispersion has 

two Dirac nodes at k = (0,0) and (7T,0), with isotropic 
dispersion in their vicinity. Including flavor indices, we 
thus have 6 two-component Dirac fermions. 

We use anti-periodic boundary conditions in both di- 
rections for L x L lattices (L even) lattices. The fc-points 
are then k x = (2n + l)n/L, n — 0...L — 1; k y = 
(2m + 1)tt/L, m = . . . L/2 — 1, which avoids zero en- 
ergy modes. Filling th negative energy states gives us 
a Slater determinant mean field wave function for each 
of the three fermion species. The Gutzwiller projected 
wave function is then easily written down as ([5]), in terms 
of this Slater determinant. Evaluating its properties 
however requires a numerical variational (determinantal) 
Monte Carlo approach] 21 ! 22 We generate a random initial 
basis state having significant overlap with the mean field 
wave function. Random pairs of sites are selected and 
updated with the Metropolis rejection rule. 10000 'ther- 
malization' sweeps(L 2 pairwise updates) are performed 
before measurements of physical quantities. Measure- 
ments are done in 100000 sweeps. The entire process is 
repeated 10 times to ensure stability of results. 

Energetics: The energy from the projected wave func- 
tion is listed in TABLE [ft] for L up to 20. We notice that 
for the 4x4 lattice our total energy — 16. 57 J is close 
to the ground state energy — 17.35J obtained in the ex- 
act diagonalization study.— More importantly, our energy 
lies below the first excited state energy —16 J obtained 
in that study, which already implies a significant overlap 
(> 42%) between our wave function and the exact ground 
state wave function. Note, our 'variational' wave function 
has no variational parameter - which makes this agree- 
ment more remarkable, especially given that there are 
24024 SU(4) singlet states already at this system size£ 

We have checked several other simple states on this 
4x4 lattice, they all show much higher energy, compared 
to the first excited state in the exact spectrum. The 
comparison is shown in TABLE IIIIl 

Wavefunction properties: We now provide evidence 
that the resulting wave function is a spin-orbital liquid. 
First, we would like to establish that it has no conven- 
tional order, to clearly show it is not a conventional state. 
Next, the specific type of liquid state being proposed - 
with emergent Dirac fermions and Z2 gauge fluxes - needs 
to be established. 
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TABLE II: Results of the projected wave function on a L x L 
square lattice with 7r-flux ansatz under anti-periodic bound- 
ary conditions in both directions. The second column shows 
energy per site. Qbox in the third column and its relation to 
box order is defined in main text. The fourth column shows 
L~ 4 scaling of spin correlation function when L is not a mul- 
tiple of four [s z (x,y) is s z at position (x, y)\. Some entries are 
empty because the numerical errors are too large. 



L 


{H)/L 2 


Qbox/L 2 


-(^(0,0)s z (L/2,0))L 4 


4 


-1.0357(4) 


2.0 


11.14(2) 


6 


-0.9238(3) 


1.7 


20.46(5) 


8 


-0.9051(2) 


1.7 


0.0(1) 


10 


-0.8995(2) 


1.6 


20.1(2) 


12 


-0.8974(2) 


1.6 


2.7(4) 


14 


-0.8966(1) 


1.6 


18.1(4) 


16 


-0.8961(1) 


1.6 




20 


-0.8956(1) 


1.6 




24 


-0.8955(1) 







TABLE III: Energy of several (variational) states for the 
SU(4) model on 4 x 4 square lattice with periodic boundary 
condition. 



state 


energy (J) 


exact ground stated 


-17.35 


exact first excited stated 


-16 


projected SO(6) Majorana fermion mean field 


-16.57 


projected SU(4) Schwinger fermion mean field 


-6.38 


orbital ferromagnetic, spin AFM state 10 ' 23 


-14.46 


box(plaquette) ordered stated - 


-12 


four-sublattice SU(4) Neel stated 






We first check for spin-orbital order. Given the SU(4) 
symmetry, it is sufficient to compute the s z correlations 
which are found to be rapidly decaying in space. The 
structure factor for various lattice size was computed, e.g. 
FIG. [5] shows the result for the L — 20 lattice. A broad 
maximum at (7r, tt) is seen, but no Bragg-peak develops. 
We thus exclude magnetic-orbital-ordering. 

A more likely order is an SU(4) singlet state that 
breaks lattice symmetry. This is the analog of the Va- 
lence Bond Solid order for SU(2) magnets. However, 
SU(4) singlets require at least four sites so box crystalline 
orders may arise. Two such natural orders were proposed 
by Li, et al.. 10 In both, the SU(4) singlets are formed on 
1/4 of the elementary plaquettes, but these are either ar- 
ranged in a square lattice, or in a body centered rectangu- 
lar lattice with aspect ratio 2. Both box orders have bond 
energies modulated at the wavevectors (tt, 0) or (tt, tt/2). 
We first check if our wave function has these correlations 
by defining E x$ = J2 P e itF (s z T z ) P ■ (s z t z ) ?+& . Then 

Obox = (E 2 p , ,) for a L x L lattice should scale as 

L 4 , if long range order is present. For example, for the 
perfect square box state which is a product state of SU(4) 
singlets, Qbox = £ 4 /36 + (13/9)L 2 . In the absence of or- 



FIG. 2: s z structure factor for the 7r-flux SO(6) projected 
wave function on a 20 x 20 square lattice with anti-periodic 
boundary conditions. 



<s z (k)s z (-k)> for L=20 




der however, this quantity will scale as L 2 . For L < 20 
we did not observe L A scaling but rather good L 2 scal- 
ing, as shown in Column 3 of TABLE [TT] indicating no 
sign or box order upto 400 site systems. An independent 
check is provided by modulating the mean field parame- 
ters Xjk to realize the box orders. The average energy of 
the projected state is then compared against the unmod- 
ulated wave function. Wc find that the energy always 
increases, for both kinds of orders, pointing to the sta- 
bility of the unmodulated state. Within mean field the- 
ory alone, however, the state is locally unstable to box 
modulation.™ The more reliable projected energy study 
however point to the opposite conclusion. 

Motivated by a recent proposal of chiral SU(N) states 
in large-N limit^ we also consider a chiral state on the 
square lattice. We add to the 7r-flux state pure imagi- 
nary hopping of fermions on the diagonal bonds in such 
a way that each triangle has +tt/2 flux. This is a par- 
ticular mass term of the Dirac fermions and it opens a 
gap in the mean field dispersion, similar to the box order 
mentioned above. Indeed the mean field energy decreases 
with the diagonal hoppings. However after projection the 
energy always increase after adding this term, indicating 
stability against this chiral order. This distinction be- 
tween mean field and projected mean field energetics has 
been observed in other projected wave function studies 
as welli^ 

We expect the spin-orbital liquid to be a nodal Z 2 
state, i.e. it contains emergent Z 2 gauge fields and nodal 
Dirac fermions that behave like free particles at low en- 
ergies. Establishing this directly is more challenging - it 
is well known that observing the Z 2 topological order of 
projected wave functions in the presence of gapless gauge 
charged fermions is tricky^ and left to future work. Free 
nodal fermions would lead to spin and orbital correla- 
tions that decay as 1/r 4 , which we check for by comput- 
ing L i (s z (0,0)s z (L/2,0)) in a size L system. The fast 
decay limits us to L < 14. The results are shown in col- 
umn 4 of TABLE [TTJ There are strong commensuration 
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FIG. 3: r z structure factor for the 7r-fiux SO(6) projected 
wave function with average fermion filling 5/8 for one fermion 
specie(other two species are still half filled), on a 20x20 square 
lattice with anti-periodic boundary conditions. 



<t z (k)x z (-k)> for L=20 




FIG. 4: (Color online) r z structure factor when k y = 7r(a ID 
cut of FIG. Hand FIG. [3]) for the tt-Aux SO (6) projected wave 
function with average fermion filling p = 1/2 (green solid line 
with symbol +) and p — 5/8(red dash line with symbol x) for 
one of the three species on a 20 x 20 square lattice with anti- 
periodic boundary conditions, both curves are of arbitrary 
scale. 
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effects which reduce the correlation when L/2 is an even 
number. However, for the other three values of L the 
correlation seems to show the required scaling. Another 
indirect evidence for such fermionization is the nature 
of these spin correlations in the presence of a Zeeman 
field AH = —hJ2 r T r> which leads to a shifted chemical 
potential for one fermion specie. We find that the pro- 
jected wave function now has a ring of incommensurate 
correlations around {-k, it) (FIG. [3] and FIG. [J). 

Breaking SU(4) symmetry It is natural to ask if 
the physical conclusions derived above are stable when 
enlarged SU(4) symmetry of our model is lost. Since the 
gauge fluxes are gapped, a weak perturbation cannot lead 
to confinement. Also, the gapless nodal fcrmions are ac- 
tually protected by discrete symmetries - one needs to 
break lattice symmetry to gap the nodes, as can be seen 
from an analysis of the fermion bilinear terms (see Ap- 
pendix [B] for details) . The only physical difference that 
arises in the lower symmetry case is that the chemical 



potential of the fermions may not be at the nodal points. 
Hence the SU(4) symmetry is not essential to our conclu- 
sions. In future^! we will apply this analysis to realistic 
Hamiltonians with reduced symmetry and search for liq- 
uid phases in those regimes. 



VI. PHYSICAL REALIZATIONS 

Since most natural spin orbit Hamiltonians are not 
near the SU(4) symmetric point, it is natural to ask where 
we might expect to find a model where the SU(4) symme- 
try is even approximately realized. As discussed in the 
introduction, cold atom systems in optical lattices pro- 
vide some promising direction for realization, if sufficient 
cooling of those magnetic Hamiltonians can be achieved. 
In this section, we point out that even in solid state sys- 
tems, approximate SU(4) symmetry may be achieved, on 
certain high symmetry lattices. In particular, we point 
out that on the diamond lattice, if only nearest neighbor 
exchange is considered, the interactions are close to the 
SU(4) point. Exchange interaction arises from a com- 
bination of hopping and on site interaction. The main 
observation is that due the high symmetry of the dia- 
mond lattice, hopping matrix elements must be SU(4) 
symmetric. The onsite interactions deviate from SU(4) 
symmetry due to, for e.g., the Hunds interaction. How- 
ever, these are typically a fraction of the overall repulsion 
leading to nearly SU(4) symmetric exchange. 

For a system of d 3z 2_ r 2 1 d x i_ y i orbitals on the dia- 
mond lattice with full lattice symmetry and without spin- 
orbital coupling, we first prove that the electron hopping 
matrix elements on nearest neighbor bonds have SU(4) 
symmetry. Denote the creation operators of the two or- 
bitals as d\ ar and d\ ar respectively, where a is spin in- 
dex, r is site index. The hopping amplitude on bond 
< ij > is generically a 2 x 2 matrix t, and the process 
is described by the term Y, a ,b,a 4<rf*ab < W Consider a 
bond from origin along the (111) direction, the reflec- 
tion x — > y, y — > x does not change this bond, however 
the orbitals transform non-trivially d a — > ^2 b (o~ z ) a bdb- 
Since this reflection is a physical symmetry, the elec- 
tronic Hamiltonian should be invariant under its action, 
thus we get t = a z ■ t ■ a z . Similarly consider a three- 
fold rotation x — ► y, y — > z, z — ► x, it does not change 
the (111) bond as well, but the orbitals transform as 
da E 6 [(-1/2)<t* + i(V3/2)o*]« 6 d6. Then we get 
t = [(-l/2)a z - i(v / 3/2)cr«] ■ t ■ [(-l/2)cr 2 + i(V3/2)a y ]. 
These two conditions on t ensure that t is proportional 
to identity matrix. Thus we have proved that the hop- 
ping on (111) direction preserves both orbital and spin, 

namely is given by the term t Y^ a ,a %ai^aaj - B ^ lattice 
symmetry we conclude that all other nearest neighbor 
bonds have this property. Therefore the nearest neighbor 
hoppings have SU(4) symmetry. One should note that 
this proof cannot be extended to next nearest neighbor 
and other generic hoppings. 
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However, the Coulomb interaction of these or- 
bitals generically does not have SU(4) symmetry. 
The onsite Coulomb interaction is given by the 
Kanamori parameters^ U J2 a n a^n a i + U' J2 a <b n a n b + 

J Ea<b,a,0 d L4p d af3 d b a + J Ea<b( d U d ll d bl d b1 + kC -)' 

and approximately U = U' + 2 J. The SU(4) symmetry 
is present only if J = and U = U 1 . We usually expect 
that J -C £/, then SU(4) is an approximate symmetry of 
the Hubbard model and thus an approximate symmetry 
of the derived spin-orbital exchange model. 

This type of systems may be realized in certain A-site 
spinels, where the magnetic ions form a diamond lattice, 
and when only the e g orbitals are active. One caveat 
is that the spinel structure allows for the next nearest 
neighbor exchange strength to be fairly large and even 
comparable to the nearest neighbor one, which may sig- 
nificantly break the SU(4) symmetry. Interestingly, the 
experimentally discussed 'spin-orbital' liquid candidate, 
FeSc2S4, is an e g system on the diamond lattice. 4 How- 
ever, it differs in two important respects from the ideal 
model considered here. First, there is a magnetic moment 
on each site that is Hunds coupled to the e g fermion, and 
second, the further neighbor exchange interactions are 
believed to be substantial in this material. 

We acknowledge support from nsf-dmr0645691 ) and 
discussions with M. Hermelc. 



APPENDIX A: PROOF OF THE EQUIVALENCE 

BETWEEN PROJECTED SO (6) MAJORANA 
MEAN FIELD STATE AND PROJECTED SU(4) 
SCHWINGER FERMION MEAN FIELD STATE 
FOR ID CHAIN. 

Consider a 4iV-site chain with periodic boundary con- 
dition. The mean field wave function for the Majorana 
fermion representation is 



2JV-1 



2N-1 



|*MF) - II c !,(2Hl)ir/(4JV) II C 2,(2fc+l)7r 



/(AN) 



k=0 
2N-1 



fe=0 



(Al) 



IT ^3,(2fc+l)7r/(4JV)|0) 
fc=0 



where |0) is fermion vacuum, c a .k is the 
Fourier transform of real space fermion op- 
erator ~c a . k (4iV)- 1 / 2 X; r Ca, r e- ifer . For 
a physically allowed real space configura- 
tion mentioned above {Uj}> { z m}) = 

n£i(4,4«) nf=i(4„A) n"=i(4 m 4 jio>, 

the overlap with the mean field wave function is 

^[{ x i\Ay 3 }A z rn}] = ({xi};{yj};{z m }\^MF) 



{z m }} ■ $[{z m }, {*<}] • $[{xi}, { yj }} 



(A2) 



where $ is the 2N x 2N Slater determinant for one 
fermion specie with the following matrix elements (p, q = 



...,2N) 



j J (2p-l)x, 

ix(2p-l)a g _' A (A3) 

e , N < q. 



Therefore the determinant is 

*[{*}, fa}] 
=(4iV)- Jv o;' ( ^ 3:i+E i« ) ]J(a; 2 * - ' 



hi 



(A4) 



i>i' 



\ I I (L0 2Xi - UJ 2X 



") YI(lu 2 v> -lu 2 ^') 

3>f 



where u> = e ,7r / ( 4Ar ) . And the overlap is 
^[{xi},{yj},{ z m}} 

^0 JJ(w 2 *™ -w 2 w) 



x ]J(lu 2 ^ 



' I I (LU 2Xi - UJ 2Zr ' 



Hi^-^n 2 n o 



u 2z ™') 2 



(A5) 



2z„ 



3>j' 



m>m' 



The mean field Hamiltonian for the standard 
'Schwinger' fermion representation is 

4 4 
#MF,a =X a a,r a a ,r+l + h ' C - ~ M YY a a,r a a,r 



r a—1 

3 



r a — 1 



r a—1 r 
3 

~ ^ mi a kr a a ,r a % a 'i-r 
r a—1 r 



(A6) 



where a' is the particle-hole conjugate of the Schwinger 
fermion a. The particle- hole transformation on the fourth 
specie is required for the following projected wave func- 
tion to represent a bosonic spin-orbital wave function. 
The quarter-filling mean field wave function for the stan- 
dard fermionic representation is (we assume N is even 
for simplicity) 

|*MF)a 

N/2-1 N/2-1 

= IT "l,(2fc+l)7r/(4JV) IT *2,(2fc+l)7r/(4JV) 
k=-N/2 k=-N/2 

N/2-1 -N/2-1 
X IT a 3 : (2fe+l)7r/(47V) II 



a' 1 

a 4,(2fc+l)7r/(4JV) 



k=-N/2 



k=-7N/2 



|0') 
(A7) 



where a a ^k is the Fourier transform of the real space 
SU(4) 'Schwinger fermions' a a ^ = (4JV) _1 / 2 e _ifcr a CT)r , 
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and |0') is the fermion 'vacuum' that can be annihilated 
by ai , a 2 , a 3 and the particle-hole conjugate of the fourth 
specie a'^. 

A physically allowed real space configuration in this 
representation is still labeled by three sets of N distinct 
numbers {#,•}, {yj}, {z m }, 

N N N 

\{xi},{vj},Um}) =Yi a \,x l Yi a iy J n Q L™ 

i—1 j—1 m—1 

N N N 



i—1 j — 1 m—1 



(A8) 



The overlap of this configuration with the mean field 
wave function is the product of four Slater determinants, 

®u[{Xi},{yj},{Zm}] = ({Xi},{yj},{ Z m}\^MF) a 
=$ a [{x t }} ■ $ a [{yi}} ■ i[{z m }] • *o'[{*iMl/i}»{*m}] 

(A9) 

The N x N Slater determinants $ a has the follow- 
ing matrix elements (p, q = 1,...,N), $ a [{xi}} pq — 

(4:N)~'-' 2 e w . The matrix elements of the 3N x 
3 N Slater determinant $ a / is (p,q = 1, . . . , 37V) 

iir(2p-l-TAf)x, 

6 4N ' q ~ N fAim 
e aw , N < q < 2N 

L e , 2N <q<3N 

Therefore the determinant of $ a is 

$ a [{xi}] = (wywj 1 -"^^ Y[(lu 2x >-^>) 

(All) 



%>%' 



And the determinant of $ a ' is 















m,i 




n (- 2% 


-^■') n -w 22 


3>j' 


m>m' 



(A12) 



field wave function is (note that ui &N = 1) 
®a[{xi},{yj},{z m }] 
=(4N)- 3N (-l) N2 u (2 - 8NH ^ x ' +zZ j 2m) 

2xi _,.,2zm\ TT7,.,2*< . .2*,/\2 



x TT(ur~" — w 



x n^-^') 2 n ( 



UJ » J — OJ JJ ) I I [CO 

j>j' m>m' 

= (-l) w V 8 ^ *»>*[{*,•}, {%•}, {z m }] 

(A13) 

Therefore these two projected wave functions for ID 
chain are identical. The crucial property utilized was 
that the Slater determinants appearing are Vandermonde 
determinants in one dimension. This property does not 
hold in higher dimensions. 



APPENDIX B: PROJECTIVE SYMMETRY 
GROUP ANALYSIS OF FERMION BILINEARS 
IN THE tt-FLUX STATE ON SQUARE LATTICE 

For the 7r-flux ansatz in FIG. [1] on an infinite lattice 
the lattice group symmetries are realized as follows (flavor 
index a omitted): 

Tx ■ C(x,y) * c (x + l,y) 

Ty : C(a; )W ) — > {-\) X C( x ^ y+ i) 

i + {-\)y - {-iy + {-iy+y 

■ {x,y) 2 y-> x ) 

Define a 4-compocnt field 

^Oi,k ^ a,ft,u' ^a.k,v'> ^a,(7r,0) + fc,v 1 ^a,(7r,0) + fc,iJ 

linearize the dispersion around the Dirac point, the low- 
energy Hamiltonian becomes 

a,|fc|«l 

where fj, are Pauli matrices acting on the u, v 2D space. 
The original fermion in real space can be represented 
as c a<f = N- it l /2 j: rkl<<1 e^{[l + (-1)»] • [ty^ + 

("Wo,*?)*] + t 1 - ' [(^ a ,g)a + (- W , S ) 3 ]}/2 

where r = (x,y). Then we have the transformation prop- 
erty of the -0 held 



T x 


^(^a; i &y ) 




Ty 


^(^a; i &y ) 




Rtt/2 


^(^a; 5 ^*-y ) 


■* (l/2)(l + ii/»)® (l+iA^H^fc.) 


m y 


r ( ^a; - ^-y ) 





Finally the overlap between this basis state and the mean 
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where v x ' v ' z are Pauli matrices acting on the 2D space 
of the two Dirac nodes (0,0) and (n, 0). The low energy 
Hamiltonian is invariant under these transformations. In 
the following we will prove that these symmetries pro- 
hibit mass term and velocity anisotropy in the low energy 
theory 

Consider a general mass term ^ Mip where M is a 4x 4 
constant non-trivial (not proportional to identity) Hermi- 
tian matrix. T x and T y translation symmetries require 
that M oc 1 ® fi with fi be a 2 x 2 constant Hermitian 
matrix. m y reflection symmetry requires n oc [i x , but this 
violate the R w /2 rotation symmetry and is forbidden. 

Consider a general velocity anisotropy term tp'(k x Mi + 
k y M2)ip with constant 4x4 matrices Mi and M 2 . Again 
T x and T y translation symmetries require that M 1( - 2 ) cx 
l<8>/ii(2). TO y reflection symmetry requires \i x \i\[i x = — H\ 
and n x H2H x — ^2, therefore /i 2 oc [i x . Also consider 180° 
rotation R\ j2 symmetry ip( k ^k y ) -» ® H v i>(-k x -k y ), 
it requires [i v [i\[i v — —(11, thus jii oc n z . Now we have 
V>t (Afo^l ® jU 2 + Bkyl ® /i 2 )^ with constants A, _B. Use 
the 90° rotation symmetry i? w / 2 we get 

A(-k y )l (g) /i z + Bfc^l <g> /i* 
=Ak x l ® [(1 - 1^)^(1 + 

+ Bfcyl ® [(1 - VK(1 + i/" y )]/2 
=Ak x l <g> /z* - Sfcyl <g> m z 

Then we have B = A and this term is just the kinetic 



energy term in the Hamiltonian. 

Consider a general bilinear term ^ M(k x , k y )ip where 
M{k Xl k y ) is a 4 x 4 hermitian matrix and a homoge- 
neous polynomial of k Xl k y . From translation symmetries 
it must be of the form 1 (8) [mo(fcE, k y )t + mi(k x , k y )fi x + 
m,2(k x , k y )[i v + m 3 (k x ,k y )fi z ] where mo, 1,2, 3 are homoge- 
neous functions of k Xl k y and are of the same order. From 
m y reflection symmetry, 

m (~k x , k y ) = m (k x , k y ), 
mi(-k x , k y ) = mi{k x , k y ), 
m 2 (-k x ,k y ) = -m 2 {k x ,k y ) 1 
m 3 (-k x ,k y ) = -m 3 (k x ,k y ) 

From i? 7r / 2 rotation symmetry, 

m (~k y , k x ) = m (k x , k y ), 
m 3 (-k y ,k x ) = -m\{k x ,k y ), 
m 2 (-k y , k x ) = m 2 (k x , k y ), 
mi{~k x , k y ) = m 3 (k x , k y ) 

This requires that too is of the form fho(k x + k y , k x k y ), 
toi is k y rhi(kl,k^), m 2 is k x k y (kl-k^)m 2 (kl + k^, k x k%), 
and to 3 is k x rhi(ky,k x ), where to ,toi,to 2 are arbitrary 
polynomials. 
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